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This note is dedicated to Professor Michele Caputo in occasion of his 10-th 
birthday. Throughout his intensive and outstanding career Professor Caputo 
has recognized the importance of the quality factor Q and fractional calculus in 
seismology, providing interesting contributions on these topics. 

Abstract 

The one-dimensional propagation of seismic waves with constant Q is shown to 
be governed by an evolution equation of fractional order in time, which interpo- 
lates the heat equation and the wave equation. The fundamental solutions for 
the Cauchy and Signalling problems are expressed in terms of entire functions 
(of Wright type) in the similarity variable and their behaviours turn out to be 
intermediate between those for the limiting cases of a perfectly viscous fluid and 
a perfectly elastic solid. In view of the small dissipation exhibited by the seis- 
mic pulses, the nearly elastic limit is considered. Furthermore, the fundamental 
solutions for the Cauchy and Signalling problems are shown to be related to 
stable probability distributions with index of stability determined by the order 
of the fractional time derivative in the evolution equation. 
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1 Introduction 



In seismology the problem of wave attenuation due to anelasticity of the Earth 
is described by the so-called quality factor Q , or, better, by its inverse Q _1 
(internal friction or loss tangent), which is related to the dissipation of the 
elastic energy during the wave propagation. Because of its great relevance in 
determining the composition and the mechanical properties of the Earth, the 
problem has been considered from different points of view by many researchers. 
Without pretending to be exhaustive, we quote (in alphabetic order of the first 
author) some original contributions and reviews among those which have at- 
tracted our attention, e.g. Aki and Richards (1980), Bcn-Mcnhahem and Singh 
(1981), Caputo (1966, 1967, 1969, 1976, 1979, 1981, 1985, 1996a) Caputo and 
Mainardi (1971), Carcione et al. (1988), Chin (1980), Futterman (1962), Gor- 
don and Nelson (1966), Jackson and Anderson (1970), Kanamori and Anderson 

(1977) , Kang and McMechan (1993) Kjartansson (1979), Knopoff (1964), Kornig 
and Muller (1989), Mitchell (1995), Murphy (1982), O'Connel and Budiansky 

(1978) , Ranalli (1987), Sabadini et al. (1985, 1987), Savage and O'Neill (1975), 
Spencer (1981), Strick (1967, 1970, 1982, 1984), Strick and Mainardi (1982), 
Yuen et al. (1986). 

It is known that seismic pulse propagation mostly occurs with a quality factor 
Q constant over a wide range of frequencies. As pointed out by Caputo and 
Mainardi (1971) and Caputo (1976), this factor turns out to be independent 
of frequency only in special linear viscoelastic media for which the stress is 
proportional to a fractional derivative of the strain, of order v less than one. 
Since these media exhibit a creep compliance depending on time by a power-law 
with exponent v, we refer to them as power-law solids, according to the notation 
by Kolsky (1956) and Pipkin (1972-1986). 

For the sake of convenience, the generalized operators of integration and 
differentiation of arbitrary order are recalled in the Appendix in the framework 
of the so-called Ricmann-Liouville Fractional Calculus. In this paper we adopt 
the Caputo definition for the fractional derivative of order a > of a causal 
function f(t) (i.e. vanishing for t < 0), 

f {m \t) if a = meN, 



f(t)-={ 1 /"* f (m) T 1-1 

dt<* J K> I — / . J . \ — dr if m - 1< a < m , V 7 

1 r(m-a) J (t-r)° ! + 1 - m 



where f( m \i) denotes the derivative of integer order m and T is the Gamma 
function. 

In Section 2 we derive the general evolution equation governing the propaga- 
tion of uniaxial stress waves, in the framework of the dynamical theory of linear 
viscoelasticity. For a power-law solid the evolution equation is shown to be of 
fractional order in time, which is intermediate between the heat equation and 
the wave equation. In fact, denoting the space and time variables by x and t 
and the response field variable by w(x,t), the evolution equation will be shown 
to be 

d 2f! w d 2 w 
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The order of the time derivative has been denoted by 2/3 for reasons that will 
appear later. Since < v < 1 , we get 1/2 < j3 < 1 . 

In Section 3 we review the analysis of the fractional evolution equation (1.2) 
in the general case < /3 < 1 , essentially based on our works, Mainardi (1994, 
1995, 1996a, 1996b). We first analyse the two basic boundary-value problems, 
referred to as the Cauchy problem and the Signalling problem, by the technique 
of the Laplace transform and we derive the transformed expressions of the re- 
spective fundamental solutions (the Green functions). Then, we carry out the 
inversion of the relevant Laplace transforms and we outline a reciprocity relation 
between the Green functions in the space-time domain. In view of this relation 
the Green functions can be expressed in terms of two interrelated auxiliary func- 
tions in the similarity variable r — \x\/( \f~Dt t3 ) . These auxiliary functions can be 
analytically continued in the whole complex plane as entire functions of Wright 
type. 

In Section 4 we show the evolution of the fundamental solutions for 1/2 < 
/3 < 1, that can be relevant in seismology to simulate the propagation of seismic 
pulses. Accounting for the low dissipation occurring in the Earth, the nearly 
elastic limit must be considered; in this case the pulse response becomes a 
narrow, sharply peaked function and the arguments by Pipkin (1972-1986) and 
Kreis and Pipkin (1986) must be adopted in order to obtain an evaluation of 
the solutions, which is suitable from numerical point of view. 

Finally, in Section 5, following Kreis and Pipkin (1986), we point out the 
interesting connection between the fundamental solution for the Signalling prob- 
lem and the density of a certain (unilateral) stable probability distribution. We 
note that this connection generalizes the one known for the standard heat equa- 
tion for which the fundamental solution for the Signalling problem is related to 
the density of the stable Levy distribution. Since the above property is expected 
to provide a further insight into our evolution equation of fractional order, the 
seismic pulse propagation with constant Q assumes an additional interest from 
a mathematical-physical point of view. 

2 Linear Viscoelastic Waves and the Fractional 
Diffusion- Wave Equation 

According to the elementary one-dimensional theory of linear viscoelasticity, the 
medium is assumed to be homogeneous (of density p) , semi- infinite or infinite in 
extent (0 < x < +oo or — oo < x < +oo) and undisturbed for t < . The basic 
equations are known to be, see e.g. Hunter (1960), Caputo & Mainardi (1971), 



Pipkin (1972-1986), Christensen (1972-1982), Chin (1980), Graffi (1982), 



Here the suffices x and t denote partial derivation with respect to space and 
time respectively, the dot ordinary time-derivation, and the star integral time- 
convolution from + to t. The following notations have been used: a for the 



o- x {x,t) = pu tt (x,t) , 
e(x,t) = u x (x,t) , 
e(x,t) = [Jo + j(t)*]a(x,t) . 



(2.1) 

(2.2) 
(2.3) 
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stress, e for the strain, J(t) for the creep compliance (the strain response to a unit 
step input of stress); the constant Jo := J(0 + ) > denotes the instantaneous 
(or glass) compliance. 

The evolution equation for the response variable w(x,t) (chosen among the 
field variables: the displacement u, the stress a, the strain e or the particle 
velocity v — Ut) can be derived through the application of the Laplace transform 
to the basic equations. We use the following notation for the Laplace transform 
of a function f(t) , locally summable for t > , 



£ {/(*)} := / ' 



f(t)dt = f(s), seC. 



and we adopt the sign to denote a Laplace transform pair, i.e. f(t) f(s) . 

We first obtain in the transform domain, the second order differential equa- 
tion 

" d 2 



in which 



dx 2 



p(s) 



w(x, s) = , 



psJ(s) 



1/2 



(2.4) 



(2.5) 



is real and positive for s real and positive. As a matter of fact, p,(s) turns out 
to be an analytic function of s over the entire s-plane cut along the negative 
real axis; the cut can be limited or unlimited in accordance with the particular 
viscoelastic model assumed. 

Wave like or diffusion like character of the evolution equation can be drawn 
from (2.5) by taking into account the asymptotic representation of the creep 
compliance for short times, 



J(t) = Jo + 0{t v ) , as t -> 0+ , 
with J > , and < v < 1 . If J > then 

hm = yfpj := - , 

s— >oo s C 



(2.6) 



(2.7) 



we have a wave like behaviour with c as the wave-front velocity; otherwise 
(Jo = 0) we have a diffusion like behaviour. In the case Jo > the wave like 
evolution equation for w(x,t) can be derived by inverting (2.4-5), using (2.6-7) 
and introducing the non dimensional rate of creep 



We get 



^■-T^>- a ' «>»■ 



/i 2 ( S ) ^s 2 [psJ(s)}^^) 2 [l + ^(s)}, 



so that the evolution equation turns out to be 

r, i / \ t d 2 W 2 2 W 

{l + ^)*}^=c 2 ^ 



(2.8) 



(2.9) 



(2.10) 
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This is a generalization of D'Alembert wave equation in that it is an integro- 
differential equation where the convolution integral can be interpreted as a per- 
turbation term. This case has been investigated by Buchen and Mainardi (1975) 
and by Mainardi and Turchetti (1975), who have provided wave-front expansions 
for the solutions. 

In the case Jo = we can re-write (2.6) as 

J(t) = ^ f ^ Y) +o(n, as^0+, (2.11) 

where, for the sake of convenience, we have introduced the positive constant 
D (whose dimensions are L 2 T v ~ 2 ) and the Gamma function T(y + 1) . Then 
we can introduce the non-dimensional function cj){t) whose Laplace transform is 
such that 

fi 2 ( S ) := s 2 [psJ{s)] = [1 + 4>(s)] . (2.12) 

Using (2.12), the Laplace inversion of (2.4-5) yields 

d 2fi w d 2 w 

where 2/3 = 2 — f so 1/2 < /3 < 1 . Here the time-derivative turns out to be 
just the fractional derivative of order 2/3 (in Caputo's sense), according to the 
Ricmann-Liouville theory of Fractional Calculus recalled in the Appendix. 

When the creep compliance satisfies the simple power-law 

J w = ^fFTir <>0 ' (2A4) 

we obtain <p(t) = 0, and the evolution equation (2.13) simply reduces to (1.2). 
As pointed out by Caputo and Mainardi (1971), the creep law (2.14) is provided 
by viscoelastic models whose stress-strain relation (2.3) can be simply expressed 
by a fractional derivative of order v . In the present notation this stress-strain 
relation reads 

<J = pD^e, Q<v<l. (2.15) 

For v = 1 the Newton law for a viscous fluid is recovered from (2.15) where D 
now represents the kinematic viscosity; in this case, since /3 = 1/2 in (1.2), the 
classical diffusion equation (or heat equation) holds for w(x, t) . When < v < 1 
the evolution equation (1.2) turns out to be intermediate between the heat 
equation and the wave equation. In general we refer to (1.2) as the fractional 
diffusion-wave equation, and its solutions can be interpreted as fractional diffu- 
sive waves, see Mainardi (1995). 

We point out that the viscoelastic models based on (2.14) or (2.15) with 
< v < 1 and henceforth governed by an evolution equation of fractional order 
in time, see (1.2) with 1/2 < j3 < 1 , are of great interest in material sciences 
and seismology. In fact, as shown by Caputo and Mainardi (1971), these models 
exhibit an internal friction independent on frequency according to the law 

Q- 1 = tan (*^pj v = - arctan (Q" 1 ) . (2.16) 
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The independence of the Q from the frequency is in fact experimentally verified 
in pulse propagation phenomena for many materials including those of seismo- 
logical interest. From (2.16) we note that the Q is also independent on the 
material constants p and D which, however, play a role in the phenomenon of 
wave dispersion. 

The limiting cases of absence of energy dissipation (the elastic energy is fully 
stored) and of absence of energy storage (the elastic energy is fully dissipated) 
are recovered from (2.16) for v = (perfectly elastic solid) and v = 1 (perfectly 
viscous fluid), respectively. 

To obtain values of seismological interest for the dissipation (Q ss 1000) we 
need to choose the parameter v sufficiently close to zero, which corresponds to a 
nearly elastic material; from (2.16) we obtain the approximate relations between 
v and Q , namely 

W0.64Q- 1 Q- 1 w 1.57 v. (2.17) 

As a matter of fact the evolution equation (1.2) turns out to be a linear 
Voltcrra integro-diffcrcntial equation of convolution type with a weakly singular 
kernel of Abel type. Equations of this kind have been treated, both with and 
without reference to the fractional calculus, by a number of authors including 
Caputo (1969, 1976, 1996b), Meshkov and Rossikhin (1970), Pipkin (1972-1986), 
Buchen and Mainardi (1975), Kreis and Pipkin (1986), Nigmatullin (1986), 
Schneider and Wyss (1989), Giona and Roman (1992), Metzler et al. (1994) 
and Mainardi (1994, 1995, 1996a, 1996b). For recent reviews on related topics 
we refer to Rossikhin and Shitikova (1997) and Mainardi (1997). 

3 The Reciprocity Relation and the Auxiliary 
Functions 




The two basic problems for our fractional wave equation (1.2) concern, for t > 0, 
the infinite interval — oo < x < +oo and the semi-infinite interval x > , 
respectively; the former is an initial - value problem, referred to as the Cauchy 
problem, the latter is an initial boundary - value problem, referred to as the 
Signalling problem. 

Extending the classical analysis to our fractional equation (1.2), and denot- 
ing by g(x) and h(t) two given, sufficiently well-behaving functions, the basic 
problems arc thus formulated as following, 

a) Cauchy problem, 

w(x, + ) = g(x) , — oo < x < +oo ; w(=foo, t) = , t>0; (3.1a) 

b) Signalling problem, 

w(x,0 + ) = 0, x>0; w(0 + ,t) = h(t) , w(+oo,t) = , t>0. (3.16) 
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If 1/2 < j3 < 1 , we must add in (3.1a) and (3.1b) the initial values of the first 
time derivative of the field variable, w t (x, + ) , since in this case (1.2) contains a 
time derivative of the second order. To ensure the continuous dependence of our 
solution with respect to the parameter (3 also in the transition from f3 = (l/2)~ 
to /? = (l/2) + , we agree to assume w t (x, + ) = . 

In view of our analysis we find it convenient from now on to add the pa- 
rameter (3 to the independent space-time variables x , t in the solutions, writing 
w = w(x, t; (3) . 

For the Cauchy and Signalling problems we introduce the so-called Green 
functions Q c {x, t; (3) and Q s {x, t; (3), which represent the respective fundamental 
solutions, obtained when g{x) — 5(x) and h(t) = S(t) . As a consequence, the 
solutions of the two basic problems are obtained by a space or time convolution 
according to 

/+oo 
g c (x-Z,t;0)g(OdS, (3.2a) 
-OO 

w{x,t;/3) = g s (x,t-T;P)h(T)dT. (3.26) 
Jo- 
lt should be noted that in (3.2a) Q c (x,t;(3) = Q c (\x\,t; (3) since the Green func- 
tion of the Cauchy problem turns out to be an even function of x. According to 
a usual convention, in (3.2b) the limits of integration are extended to take into 
account for the possibility of impulse functions centred at the extremes. 

For the standard diffusion equation (/? = 1/2) it is well known that 

g c (x, i; 1/2) := g d c (x, t) = — L= t" V=» e^ 2 /(4 D t) ^ 

g s (x,t-,i/2) ~gf(x,t) = ^^ t -m c -x 2 /(±Dt) (336) 

2v nD 

In the limiting case {3 = 1 we recover the standard wave equation, for which, 
putting c = \[T) , 

g c {x, t; 1) := g?(x, t) = ^ [6{x - ct) + 6(x + ct)] , (3.4a) 

g s (x,t;l) :=g™{x,t) = 5{t-x/c). (3.46) 

In the general case < (3 < 1 the two Green functions will be determined by 
using the technique of the Laplace transform. This technique allows us to obtain 
the transformed functions g c (x, s; (3), g s (x, s; (3), by solving ordinary differential 
equations of the 2-nd order in x and then, by inversion, g c (x, t; (3) and g s (x, t; (3). 

For the Cauchy problem (3.1a) the application of the Laplace transform 
to (1.2) with w(x,t) = g c (x,t;(3) leads to the non homogeneous differential 
equation satisfied by the image of the Green function, g c (x, s; (3) , 

D d ^-s 2l3 g' c = -6(x)s 20 - 1 , -oo<:r<+oo. (3.5) 



7 



Because of the singular term 6(x) we have to consider the above equation sepa- 
rately in the two intervals x < and x > 0, imposing the boundary conditions at 
x = =Foo , Gc{T°°, t;/3) = , and the necessary matching conditions at x = ± . 
We obtain 

cT c (x,s;(3) = .} C -(N/V^K, -oo<.t<+oo. (3.6) 



For the Signalling problem (3.1b) the application of the Laplace transform 
to (1.2) with w(x, t) = Q s (x, t; (3) leads to the homogeneous differential equation 

D^-^Gl^O, x>0. (3.7) 

Imposing the boundary conditions at x = , G s (0 + 7 t;(3) = h(t) = 5(t) , and at 
x = +oo , Q s (+oo, t; ft) =0, we obtain 

g s (x, S ;P)=c-( x /^ D ) s \ x>0. (3.8) 

From (3.6) and (3.8) we recognize for the original Green functions the following 

reciprocity relation 

2/3xG c (x,t;f3) =tg s (x,t;(3) , x>0, t>0. (3.9) 

This relation can be easily verified in the case of standard diffusion (/? = 1/2), 
where the explicit expressions (3.4) of the Green functions leads to the identity 

xG d c (x,t) =tg d s {x,t) = ^ r ^=- c- x2 /^ Dt ) = F d (r) = - M\r) , (3.10) 
2V 71 " VDt 2 

where r — xj{\J~D t 1 / 2 ) > is the well-known similarity variable and 

M d (r) = ^e- r2 /4. (3.11) 
v 7r 

We refer to F d (r) and M d (r) as to the auxiliary functions for the diffusion 
equation because each of them provides the fundamental solution through (3.10). 
We note that M d (r) satisfies the normalization condition J^°M d (r) dr = 1. 

Applying in the reciprocity relation (3.9) the complex inversion formula for 
the transformed Green functions (3.6) and (3.8), and changing the integration 
variable in a = s t , we obtain 

2/3 x g c (x, t;p)=t G s {x, t- /3) = F(r; (3) = [3r M (r; /?) . (3.12) 

where _ 

r = x/{^fDt' 3 ) >0 (3.13) 

is the similarity variable and 

: = h 1/ " " ■ =- h if - <««> 

are the auxiliary functions. In (3.14) Br denotes the Bromwich path and r > , 
< [3 < 1. 



8 



The above definitions of F(r; f3) and M(r; j3) by the Bromwich representation 
can be analytically continued from r > to any z € C, by deforming the 
Bromwich path Br into the Hankel path Ha , a contour that begins at a = 
— oo — ia (a > 0), encircles the branch cut that lies along the negative real axis, 
and ends up at a = — oo + ib (b > 0). 

The integral and series representations of F(z;8) and M(z;/3), valid on all 
of C , with < 8 < 1 turn out to be 



F(z;8) 



= E 



M(z;B) 



2™ J Ha 

{- z y 



— za^ 



da 



^ nlT(-dn) 

n=l ^ ' 



(3.15) 



E 

n=l 



-z) r 



ra! 



r(^n + 1) sin(7r/3n) 



£ do- 



71=0 



n!r[-/3n+(l-/3)] 



(3.16) 



i (-*)' 

7T ^ (n- 

n— 1 v 



1)! 



T(/?n) sin(7r/?n) 



In the theory of special functions, see Erdelyi (1955), we find an entire 
function, referred to as the Wright function, which reads (in our notation) 



2 ™ J Ha 



erf 



n -j„!r(An + /i) ' 



ze C. 



(3.17) 



where A > — 1 and \i > 0. From a comparison among (3.15-16) and (3.17) 
we recognize that the auxiliary functions are related to the Wright function 
according to 

F(z;8) = W-p fi (-z) = 8zM(z;8), M(z; 8) = W-p,i- P {-z) . (3.18) 

Although convergent in all of C, the series representations in (3.15-16) can 
be used to provide a numerical evaluation of our auxiliary functions only for 
relatively small values of r , so that asymptotic evaluations as r ^ +oo are 
required. Choosing as a variable r/8 rather than r , the computation by the 
saddle-point method for the M function is easier and yields, see Mainardi and 
Tomirotti (1995), 



r (/?-l/2)/(l-0) 
M(r/8; 8) = exp 
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AM 1 -?) 



r^+oo. (3.19) 



We note that the saddle-point method for 8 = 1/2 provides the exact result 
(3.11), i.e. M(r; 1/2) = M d (r) = cxp{-r 2 /4)/yfir , but breaks down for 8 -> 1". 
The case 8 = 1 , for which (1.2) reduces to the standard wave equation, is of 
course a singular limit also for the series representation since M(r; 1) = S(r — 1). 
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The exponential decay for r — > +00 ensures that all the moments of M(r; 0) 
in R + arc finite; in particular, see Mainardi (1997), we obtain 

r + °° r(n + i) 

/ r n M(r;0)dr = \ n = 0,l,2... (3.20) 

Jo r (P« + 1) 

In Fig. 1 we exhibit plots of the auxiliary function M(r; 0) in < r < 4 for some 
rational values of . The plots are obtained by means of a numerical matching 
between the series and the saddle-point representations. As a matter of fact it 
turns out that the function M(r;0) is decreasing for < < 1/2, while for 
l/2</3<lit first increases and then decreases exhibiting the maximum value 
M o (0) at a certain point r n (/3); as — > 1~ , Mo(0) — > +00 and r n (/3) — > 1 . 




Figure 1: Comparison of M(r; /?) (continuous line) with M(r; 1/2) (dashed line) 
in < r < 4, for various value of /?; (a) 1/4, (b) 1/3, (c) (2/3), (d) 3/4. 



4 The Evolution of the Seismic Pulse from the 
Fundamental Solutions 

It is known that in theoretical seismology the delta-Dirac function is of great 
relevance in simulating the pulse generated by an ideal seismic source, concen- 
trated in space (S(x)) or in time (5(t)). Consequently, the fundamental solutions 
of the Cauchy and Signalling problems are those of greater interest because they 
provide us with information on the possible evolution of the seismic pulses dur- 
ing their propagation from the seismic source. Accounting of the reciprocity 
relation (3.12) and the similarity variable (3.13), r = x/(\/~Dt^) , the two fun- 
damental solutions can be written, for x > and t > , as 

Gc(x, t; p) = ^ F(r; p) = M(r; (3) , (4.1a) 

Gs(x,t;f3) = if (r;/?) = J* M(r; 0) . (4.16) 
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The above equations mean that for the fundamental solution of the Cauchy 
[Signalling] problem the time [spatial] shape is the same at each position [in- 
stant] , the only changes being due to space [time] - dependent changes of width 
and amplitude. The maximum amplitude in time [space] varies precisely as 1/x 
[1/t] . The two fundamental solutions exhibit scaling properties that make easier 
their plots versus distance (at fixed instant) and versus time (at fixed position). 
In fact, using the well-known scaling properties of the Laplace transform in (3.6) 
and (3.8), we easily prove, for any p . q > , that 

g c (p Xl qt; P) = ~ Qc(px/q ,t; 0) , (4.2a) 
qP 

Q s {px,qt-P) = -Q s {px/qf , ,t;p), (4.26) 

and, consequently, in plotting we can choose suitable values for the fixed vari- 
able. 




Figure 2: Comparison of the representations of M (r; /3) with [3 = 1—e around the 
maximum r m 1 , in the cases (a) e = 0.01 , (b) e = 0.001 , obtained by Pipkin's 
method (continuous line), 100 terms-series (dashed line) and and saddle-point 
method (dashed-dotted line). 



In order to inspect the evolution of the initial pulse for seismological pur- 
poses, we need to plot Q c (x,t;P) versus x and Q s {x,t;0) versus t as f3 is suffi- 
ciently close to 1 (nearly elastic cases) to ensure a sufficiently low value for the 
constant internal friction Q _1 . From (2.13)-(2.14) and (2.16)-(2.17) we need to 
consider /3 = 1 — e with e = i//2 of the order of 0.001 to 0.01 . In the evaluation of 
the auxiliary functions in the nearly elastic cases, we note that the matching be- 
tween the series and saddle point representations is no longer achieved since the 
saddle point turns out to be wide and the consequent approximation becomes 
poor. In these cases we need to adopt the ingenious variant of the saddle-point 
method introduced by Pipkin (1972-1986), see also Kreiss and Pipkin (1986), 
which allows us to see some structure in the peak while it tends to the Dirac 
delta function. With Pipkin's method we get the desired matching with the 
series representation just in the region around the maximum r ~ 1 , as shown 
in Fig. 2a, b, where we exhibit the significant plots of the auxiliary function 
M(r; (i) with (3 = 1 - e for e = 0.01 and e = 0.001 . 
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Figure 4: Plots of the fundamental solution G s (x, t; (3) versus t at fixed x = D = 
1 , with /3 = 1 - e in the cases (a) e = 0.01 , (b) e = 0.001 . 



Once obtained the auxiliary function M(r; (3) in the nearly elastic cases, we 
easily get the corresponding plots of the fundamental solutions of the Cauchy 
and Signalling problems by using (4.1a-b), see Figs. 3a,b and 4a, b. 

We also note the exponential decay of Q c (x,t;(3) as x — > +oo (at fixed t) 
and the algebraic decay of G s ( x > *j /3) as t — > +oo (at fixed x), for < j3 < 1 . In 
fact, using (4.1a,b) with (3.16) and (3.19), we get 



exp 



-b(t)x 



1/(1-/9) 



00 , 



G s (x,t;P) - c{x)t- {1+f3) , t^oo. 
where a(t) , b(t) and c{x) are positive functions. 



(4.3a) 
(4.36) 
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5 The Fundamental Solutions as Probability 
Density Functions 

It is well known that the fundamental solution of the standard diffusion equation 
for the Cauchy problem is related with the Gauss or normal probability law, 
bilateral in space. In fact, recalling (3.3a), we have 



G d c (x,t) = 



1 



2VirDt 



where 



p G (x;a) :=-^ e -- 2 /(^ 2 ), (7 ^2Dt, 



(5.1) 



(5.2) 



denotes the well-known Gauss probability density function (pdf) with variance 
a 2 . The associated cumulative distribution function (calf ) is known to be 



x;a) := J 



p G (x;a)dx 



1 + erf 



'2(7 



1 + erf 



2 VDt 



(5.3) 



The moments of even order of the Gauss pdf turn out to be with neN, 



+°° (2nV 

x 2n p G (x; cr) dx = a 2n = (2n - 1)!! a 2n = (2n - 1)!! (2Dt) n . (5.4) 



If we consider the fundamental solution of the standard diffusion equation 
but for the Signalling problem, we note that it is related to the Levy probability 
law, unilateral in time (a property not so well-known as that for the Cauchy 
problem!). In fact, recalling (3.3b), we have 



2\Ar^Di 3 / 2 



-x 2 /(4Dt) 



where 



PL(t;fJ.) = 



VJ 1 



-/V(2t) u = 



/2^t 3 / 2 

denotes the Levy pdf , see Feller (1971), with cdf 



2D ' 



Jo 



dt = erfc 



crfc 



2 VDt 



(5.5) 



(5.6) 



(5.7) 



The Levy pdf has all moments of integer order infinite, since it decays at infinity 
as i -3 / 2 . However, we note that the moments of real order 5 are finite only if 
0<(5<l/2.In particular, for this pdf the mean (expectation) is infinite, 
but the mediane is finite. In fact, from V L{tmed\ I*) = 1/2, it turns out that 
tmed ~ 2/i, since the complementary error function gets the value 1/2 as its 
argument is approximatively 1/2 (a better evaluation of the argument is 1/2.1). 

The Gauss and Levy laws are special cases of the important class of a - 
stable probability distributions, or stable distributions with index of stability 
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(or characteristic exponent) a = 2 and a = 1/2, respectively. Another special 
case is provided by the Cauchy law with pdf pc(x;X) = X/[n(x 2 + A 2 )] and 
a = 1. 

The name stable has been assigned to these distributions because of the 
following property: if two independent real random variables with the same 
shape or type of distribution arc combined linearly and the distribution of the 
resulting random variable has the same shape, the common distribution (or its 
type, more precisely) is said to be stable. More precisely, if Y\ and Y 2 are random 
variables having such distribution, then Y defined by the linear combination 
cY = ci Yi + C2 Y-i has a similar distribution with the same index a for any 
positive real values of the constants c , C\ and c 2 with c a = c" + c% . As a matter 
of fact only the range < a < 2 is allowed for the index of stability. The case 
a = 2 is noteworthy since it corresponds to the normal distribution, which is 
the only stable distribution which has finite variance, indeed finite moments of 
any order. In the cases < a < 2 the corresponding pdf p a (y) have inverse 
power tails, i.e. Jj . >A p a (2/) dy = 0(\~ a ) and therefore their absolute moments 
of order S are finite if < S < a and infinite if S > a . 

The inspiration for systematic research on stable distributions, originated 
with Paul Levy, was the desire to generalize the celebrated Central Limit The- 
orem (CLT). 

The restrictive condition of stability enabled some authors to derive the 
general form for the characteristic function (cf, the Fourier transform of the 
pdf ) of a stable distribution, see Feller (1971). A stable cf is also infinitely 
divisible, i.e. for every positive integer n it can be expressed as the nth power 
of some cf. Equivalently we can say that for every positive integer n a stable 
pdf can be expressed as the n-fold convolution of some pdf . All stable pdf are 
unimodal and indeed bell-shaped, i.e. their n-th derivative has exactly n zeros, 

Using standardized random variables, the a-stable distributions turn out to 
depend on an additional parameter 7 , said the skewness parameter. Denoting 
a stable pdf by p a (y; 9) , we note the property p a (—y; —0) = p a (y; 9) . Conse- 
quently a stable pdf with 6 = is necessarily symmetrical. As a matter of fact 
|0| < a if < a < 1 and \9\ < 2 - a if 1 < a < 2 . Stable distributions with 
extremal values of 9 are called extremal. 

From the theory one recognizes that the normal distribution is the only stable 
df independent on 9, and that all the extremal stable distributions with < 
a < 1 are unilateral, i.e. vanishing in IR* if 9 = ±a . In particular, the following 
representations by convergent power series are valid for stable distributions with 
< a < 1 (negative powers) and 1 < a < 2 (positive powers) , for y > , 



Pa(v, 0) = — f2(-y~ a T F(na , + 1) sin '^(9 -a) , < a < 1 , (5.8) 
7r y n! 12 J 

y 71=1 

Pa (y,9) = — f^(-yr F(n/a + 1) S in&g-a)l , Ka<2. (5.9) 



niT 



From (5.8)-(5.9) a relation between stable pdf with index a and 1/a can be 
derived. Assuming 1/2 < a < 1 and y > , we obtain 

yjSt+i Pi/«(jT a ; ) = Pa(y; 0*) , 0* = a(9 + 1) - 1 . (5.10) 
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A quick check shows that 9* falls within the prescribed range, \0*\ < a , provided 
that \9\ < 2 — 1/a. Furthermore, we can derive a relation between extremal 
stable pdf and our auxiliary functions of Wright type. In fact, by comparing 
(5.8-9) with the series representations in (3.15-16) and using (3.18), we obtain 

1 a 

p a (y; -a) = - F(y- a :a) = — - M(y- a ;a) , < a < 1 , (5.11) 

y y a+L 

Pa (y; a - 2) = - F(y; 1/a) = - M(y; 1/a) ,Ka<2. (5.12) 

y & 

Consequently we can interpret the fundamental solutions (4.1a) and (4.1b) in 
terms of stable pdf, so generalizing the arguments for the standard diffusion 
equation based on (5.1-7). 

We easily recognize that for < (3 < 1 the fundamental solution for the 
Signalling problem provides a unilateral extremal stable pdf in (scaled) time 
with index of stability a = f3 , which decays according to (4.3b) with a power 
law. In fact, from (4.1b) and (5.11) we note that, putting y = r -1 /' 3 = t , 

(x/y/D) 1 /Pg s (x,t;p)=p (T;-(3), r = t (Vd/x) 1 " 3 > . (5.13) 

This property has been noted also by Kreiss and Pipkin (1986) based on (3.8) 
and on Feller's result, p a (t; —a) 4- exp(— s a ) for < a < 1 . 

As far as the Cauchy problem is concerned, we note that the corresponding 
fundamental solution provides a bilateral symmetrical pdf in (scaled) distance 
with two branches, for x > and x < , obtained one from the other by 
reflection. For large |x| each branch exhibits an exponential decay according to 
(4.3) and, only for 1/2 < (3 < 1 , it is the corresponding branch of an extremal 
stable pdf with index of stability a = 1/(3 . In fact, from (4.1b) and (5.12) we 
note that, putting y = \r\ = £ > , 

20VDtf'g c (\x\,t;l3)=p 1/ f ) (S,l/l3-2), £ = \x\/(VDtP) > . (5.14) 

This property had to the authors' knowledge not been noted: it properly gener- 
alizes the Gaussian property of the pdf found for f3 = 1/2 (standard diffusion). 
Furthermore, using (3.20), the moments (of even order) of Q c (x,t;j3) turn out 
to be 

£°°x 2n G c (x, t; 13) dx - r(2 2 (3n + +i) ( ^ } " ' n = ° ' 1 ' 2 ' ' ' ' (5A5) 

We recognize that the variance is now proportional to Dt 2 @ , which implies a 
phenomenon of fast diffusion ifl/2</3<l. 

Appendix: Essentials of Fractional Calculus 

Fractional calculus is the field of mathematical analysis which deals with the 
investigation and applications of integrals and derivatives of arbitrary order. 
The term fractional is a misnomer, but it is retained following the prevailing 
use. 
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According to the Ricmann-Liouville approach to fractional calculus, the no- 
tion of fractional integral of order a (a > 0) is a natural consequence of the well 
known formula (usually attributed to Cauchy), that reduces the calculation of 
the n— fold primitive of a function f(t) to a single integral of convolution type. 
In our notation the Cauchy formula reads 

J n f{t)--=Ut) = T ^— /"(t-r)"- 1 f(T)dr, t>0, neN, (Al) 

where IN is the set of positive integers. From this definition we note that f n (t) 
vanishes at t = with its derivatives of order 1, 2, . . . , n — 1 . For convention 
we require that f(t) and henceforth f n (t) be a causal function, i.e. identically 
vanishing for t < . 

In a natural way one is led to extend the above formula from positive integer 
values of the index to any positive real values by using the Gamma function. 
Indeed, noting that (n — 1)! = T(n) , and introducing the arbitrary positive real 
number a , one defines the Fractional Integral of order a > : 

J«/(t):=-L [ (t-rr- 1 f(r)dr, t>0, a € R + , (A2) 
T(a) Jo 

where H + is the set of positive real numbers. For complementation we define 
J° := / (Identity operator), i.e. we mean J° f(t) = f(t). Furthermore, by 
J a f{0 + ) we mean the limit (if it exists) of J a f(t) for t — > + ; this limit may 
be infinite. 

We note the semigroup property J a J^ = J a+ ^ , a, (3 > , which implies 
the commutative property J@J a — J a J^ 7 and the effect of our operators J a on 
the power functions 

J a f = = F(7 1 ~ x tl+a > «>0, 7>-l. *>0. (A3) 

r( 7 + i + a) v ; 

These properties are of course a natural generalization of those known when the 
order is a positive integer. 

Introducing the Laplace transform by the notation C {/(£)} := 
/ °°e _s * f(t) dt = f(s) , s e C , and using the sign to denote a Laplace trans- 
form pair, i.e. f(t) f(s) , we note the following rule for the Laplace transform 
of the fractional integral, 

J a f(t) + ^, a>0, (A4) 

which is the generalization of the case with an n-fold repeated integral. 

After the notion of fractional integral, that of fractional derivative of order 
a (a > 0) becomes a natural requirement and one is attempted to substitute a 
with —a in the above formulas. However, this generalization needs some care in 
order to guarantee the convergence of the integrals and preserve the well known 
properties of the ordinary derivative of integer order. 

Denoting by D n with n € N , the operator of the derivative of order n , we 
first note that D n J" = I , J™ D n ^ I , n e IN , i.e. D n is left-inverse (and 
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not right-inverse) to the corresponding integral operator J" . In fact we easily 
recognize from (A.l) that 



71-1 



J"£>"/(i) = /(t)-^/( fe )(0+)-, t>0. 



k=0 



(A5) 



As a consequence we expect that D a is defined as left-inverse to J a . For this 
purpose, introducing the positive integer to such that to — 1 < a < to , one 
defines the Fractional Derivative of order a > as D° f(t) := D m J m ^ a f(t) ; 
i.e. 



d n 
dt r ' 



1 



f(r) 



T(m-a) Jo (t-r)^ 1 -™ 
/(*), 



m — 1 < a < to, 



(A6) 



a = m . 



Defining for complementation D° = J° = I , then we easily recognize that 
D a J a = I, a > , and 



D a t 1 



r(7 + i) 



a > 0, 7 > -1, t > 0. 



(A.7) 



r(7 + 1 - a) 

Of course, these properties are a natural generalization of those known when 
the order is a positive integer. 

Note the remarkable fact that the fractional derivative D a f is not zero for 
the constant function f(t) = 1 if a £ IN . In fact, (A. 7) with 7 = teaches us 
that 

t~ a 



D a l = 



a > 0, t > 0. 



(A.8) 



r(l-a) 

This, of course, is = for a € IN, due to the poles of the gamma function in the 
points 0, —1, —2, . . .. We now observe that an alternative definition of fractional 
derivative, originally introduced by Caputo (1967) (1969) in the late sixties and 
adopted by Caputo and Mainardi (1971) in the framework of the theory of Linear 
Viscoelasticity, is D" f(t) := J m ~ Q D m f(t) with to - 1 < a < m , melN, i.e. 

rt /(m)( T ) 



D?f(t) := 



1 

T(m-a) J (t - r) a + l - 
d m 

T-/(*)- 
dt m 



dr , m — 1 < o; < to , 



(A9) 



a = to . 



This definition is of course more restrictive than (A. 6), in that requires the 
absolute integrability of the derivative of order to. Whenever we use the operator 
D" we (tacitly) assume that this condition is met. We easily recognize that in 
general 



D a f(t) := D m J m - a f(t) ^ J m ~ a D m f(t) := D* f{t) 



(A10) 



unless the function f(t) along with its first to — 1 derivatives vanishes at t = + . 
In fact, assuming that the passage of the m-derivative under the integral is 
legitimate, one recognizes that, for to — 1 < a < to and t > , 



m— 1 



D a f(t) = D«f(t)+J2 



k=0 



T(k-a + l) 



-/W(0+), (All) 
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and therefore, recalling the fractional derivative of the power functions (A. 7), 

/ rn-l k \ 

D a f /(t) - £ - /< fc >(0+)J = /(t) . (A12) 

The alternative definition (A. 9) for the fractional derivative thus incorporates 
the initial values of the function and of its integer derivatives of lower order. The 
subtraction of the Taylor polynomial of degree to — 1 at t = + from f(t) means 
a sort of regularization of the fractional derivative. In particular, according to 
this definition, the relevant property for which the fractional derivative of a 
constant is still zero can be easily recognized, i.e. 

D°1 = 0, a>0. (A13) 



We now explore the most relevant differences between the two fractional 
derivatives (A. 6) and (A. 9). We agree to denote (A. 9) as the Caputo frac- 
tional derivative to distinguish it from the standard Riemann-Liouville frac- 
tional derivative (A. 6). We observe, again by looking at (A. 7), that D a t a ~ x = 
, a > , t > . From above we thus recognize the following statements about 
functions which for t > admit the same fractional derivative of order a , with 
to — 1 < a < m , to G IN , 

m 

D a f(t) = D a g(t) ^ f(t) = g{t) + ]T c 3 t a ~i , (A14) 

m 

0? f(t) = D" git) f{t) = git) + ]T Cj t m ~i . (A15) 

In these formulas the coefficients Cj are arbitrary constants. 

For the two definitions we also note a difference with respect to the formal 
limit as a — > (to — 1) + . From (A. 6) and (A. 9) we obtain respectively, 

a — )• (to — 1) + D a fit) -^D m J fit) = D™- 1 fit) ; (A16) 

a -+ (to - 1)+ => /(t) J D m /(t) - D™- 1 fit) - Z^" 1 ) (0+) . (A17) 



We now consider the Laplace transform of the two fractional derivatives. 
For the standard fractional derivative D a the Laplace transform, assumed to 
exist, requires the knowledge of the (bounded) initial values of the fractional 
integral J m ~ Q and of its integer derivatives of order k = 1, 2, . . . , to — 1 . The 
corresponding rule reads, in our notation, 

m— 1 

D a fit) + s a fis) -J2 Dk j(m ~ a) ^ 0+ ) s m_1_fe , m - 1< a < m . (A18) 

fc=0 

The Caputo fractional derivative appears more suitable to be treated by 
the Laplace transform technique in that it requires the knowledge of the 
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(bounded) initial values of the function and of its integer derivatives of order 
k = 1,2,..., tyl — 1 , in analogy with the case when a = m . In fact, by using 
(A. 4) and noting that 

m—1 j, 

J a D? f(t) = J a J m - a D m f(t) = J m D m f(t) = f(t) - ]T / (fe) (0 + ) M - (A19) 
we easily prove the following rule for the Laplace transform, 

m — 1 

D« f{t)-rs a J(s)- ]T /W(0+) S a - 1 - fe , m-l<a<m. (A.20) 

Indeed, the result (A.20), first stated by Caputo (1969) by using the Fubini- 
Tonelli theorem, appears as the most "natural" generalization of the corre- 
sponding result well known for a — m . 

This appendix is based on the review by Gorenflo and Mainardi (1997). For 
more details on the classical treatment of fractional calculus the reader is referred 
to Erdelyi (1954), Oldham and Spanier (1974), Samko et al. (1987-1993) and 
Miller and Ross (1993). Gorenflo and Mainardi have pointed out the major util- 
ity of the Caputo fractional derivative in the treatment of differential equations 
of fractional order for physical applications. In fact, in physical problems, the 
initial conditions are usually expressed in terms of a given number of bounded 
values assumed by the field variable and its derivatives of integer order, no 
matter if the governing evolution equation may be a generic integro-differential 
equation and therefore, in particular, a fractional differential equation. 
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